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Abstract — In this paper, we investigate the Hermite spectral 
method (HSM) to numerically solve the forward Kolmogorov 
equation (FKE). A useful guideline of choosing the scaling 
factor of the generalized Hermite functions is given in this 
paper. It greatly improves the resolution of HSM. The conver- 
gence rate of HSM to FKE is analyzed in the suitable function 
space and has been verified by the numerical simulation. 
As an important application and our primary motivation to 
study the HSM to FKE, we work on the implementation 
of the nonlinear filtering (NLF) problem with a real-time 
algorithm developed in 1 17 1. The HSM to FKE is served as 
the off-line computation in this algorithm. The translating 
factor of the generalized Hermite functions and the moving- 
window technique are introduced to deal with the drifting 
of the posterior conditional density function of the states in 
the on-line experiments. Two numerical experiments of NLF 
problems are carried out to illustrate the feasibility of our 
algorithm. Moreover, our algorithm surpasses the particle 
filter as a real-time solver to NLF. 



I. Introduction 

The central problem in the field of nonlinear filtering 
(NLF) is to give the instantaneous and accurate estimation 
of the states based on the noisy observations, if enough 
computational resources are provided. In 1960s, Duncan 
fTl . Mortensen (T8| and Zakai |23| independently derived 
the so-called DMZ equation, which the unnormalized 
conditional density function of the states satisfies. Hence 
the central problem in NLF is translated into solving the 
DMZ equation in the real time and memoryless manner. It 
is worthy to point out that the "real time" and "memory- 
less" are the most important properties one would like to 
maintain in the design of the optimal/suboptimal nonlinear 
filters for real applications. More specifically, "memory- 
less" refers that one only needs the latest observation to 
update the estimation of the states without refering back 
to any earlier observation history; "real time" means that 
the decision of the states is made on the spot, while the 
observation data keep coming in. 

It is well known that the exact solution to the DMZ 
equation, generally speaking, can not be written in a 
closed form. With the well-posedness theory of the DMZ 
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equation in mind, many mathematicians inake efforts to 
seek an efficient algorithm to construct a "good" approxi- 
mate solution to the DMZ equation. One of the methods is 
the splitting-up method from the Trotter product formula, 
which was first described in Besoussan, Glowinski, and 
Rascanu [5], |6|. It has been extensively studied in many 
articles later, for instance 1121 . 1141 . 1151 and |19| . In 
1990s, Lototsky, Mikulevicius and Rozovskii fT6) devel- 
oped a new algorithm (so-called S'^-algorithm) based on 
the Cameron-Martin version of Wiener chaos expansion. 
However, the above methods require the boundedness of 
the drifting term and the observation term (/ and h in 
l |l.l[ )), which leaves out even the linear case. To overcome 
this restriction, Yau and Yau |22| developed a novel 
algorithm to solve the "pathwise-robust" DMZ equation, 
where the boundedness of the drift term and observation 
term is replaced by some mild growth conditions on / 
and h. Nevertheless, they still made the assumption that 
the drift term, the observation term and the diffusion term 
are "time-invariant". That is to say, /, h and g in ([TTTJ 
are not explicitly time-dependent. In 1171 . we generalized 
Yau-Yau's algorithm to the most general settings of the 
NLF problems - "time-varying" case - where /, h and g 
are explicitly time-dependent. 

Our study of solving the forward Kolmogorov equation 
(FKE) by the Hermite spectral method (HSM) is closely 
related to the implementation of the algorithm developed 
in 1171 . The detailed formulation of our algorithm could 
be found in Appendix A or 1 17|. Briefly speaking, in our 
algorithm, we start from the signal based model: 

dxt = f{xt,t)dt + g{xt, t)dvt, 
dyt = h{xt,t)dt + dwt, 

where Xt is a vector of the states of the system at time 
t with xo satisfying some initial distribution and yt is 
a vector of the observations at time t with yo = 0. 
vt and wt are vector Brownian motion processes with 

E[dvtdvf] = Q{t)dt and E[dwtdw'^] = S{t)dt, S{t) > 
0, respectively. The DMZ equation is derived as 

da(x, t) = La{x, t)dt + a[x, t)h^ {x, t)S^^ {t)dyt 
a{x, 0) = ao{x), 

(1.2) 

where a{x) is the unnormalized conditional density funci- 
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ton, (Jo{x) is the density of the initial states xq, and 



dxidxj 



dxi 



(1.3) 



To maintain the real-time property, solving the DMZ 
equation is translated into solving a FKE off-line and 
updating the initial data on-line at the beginning of each 
time interval. Let Vk = {0 = ro < ri < ■ ■ ■ < ^ T} 
be a partition of [0,r]. The FKE needs to be solved at 
each time step is 

dui , 

(1.4) 

where L is defined as l |1.3[ l. The initial data is updated as 
follows 



-(x, t) = - ^/i^S ^h^Ui{x,t) on[Ti-i,Ti], 



ui{x,0) = ao{x) 



Wi(a;,Ti_i) = exp ri_i)5 ^(ri_i) 



i > 2, 



_ilx,ri_il 



(1.5) 



where Ui is transformed from o, see the detailed for- 
mulation of our algorithm in the Appendix A or 1171 . 
From the above description, it is not hard to see that the 
FKE \\A\ needs to be solved repeatedly on each time 
interval [ri_i,ri]. Thus, it is crucial to obtain a good 
approximate solution to \\A) . In this paper, we adopt 
HSM to solve FKE for two reasons: on the one hand, 
HSM is particularly suitable for functions defined on the 
unbounded domain which decays exponentially at infinity; 
on the other hand, HSM could be easily patched with 
the numerical solution obtained in the previous time step 
while the moving-window technique is in use in the on- 
line experiments. 

The HSM itself is also a field of research, which 
could be traced back to 1970s. In |10|, Gottlieb et. al. 
gave the example sin a; to illustrate the poor resolution 
of Hermite polynomials. To resolve M wavelength of 
sin a;, it requires nearly K'P Hermite polynomials. Due 
to this fact, they doubted the usefulness of Hermite 
polynomials as bases. The Hermite functions inherit the 
same deficiency from the polynomials. Moreover, it is lack 
of fast Hermite transform (some analogue of fast Fourier 
transform). Despite of all these drawbacks, the HSM has 
its inherent strength. Many physical models need to solve 
a differential equation on an unbounded domain, and 
the solution decays exponentially at infinity. From the 
computational point of view, it is hard to describe the 
rate of decay at infinity numerically or to impose some 
artificial boundary condition cleverly on some faraway 
"boundary". Therefore the Chebyshev or Fourier spectral 
methods are not so useful in this situation. As to the HSM, 
how to deal with the behavior at infinity is not necessary. 
Recent applications of the HSM can be found in [S], (51 , 
fTTl . (20J, tlJJ, etc. 



To overcome the poor resolution, a scaling factor is 
necessary to be introduced into the Hermite functions, 
refer to fS), (4). It is shown in [4| that the scaling factor 
should be chosen according to the truncated modes A'^ 
and the asymptotical behavior of the function f{x), as 
\x\ — >■ oo. Some efforts have been made in seeking 
the suitable scaling factor a, see 

II, QU, (20), etc. To 
optimize the scaling factor is still an open problem, even 
in the case that fix) is given explicitly, to say nothing 
of the exact solution to a differential equation, which is 
generally unknown a-priori. Although some investigations 
about the scaling factor have been made theoretically, 
as far as we know, there is no practical guidelines of 
choosing a suitable scaling factor. Nearly all the scaling 
factors in the papers with the application of HSM are 
obtained by the trial-and-error method. Thus, we believe 
it is necessary and useful to give a practical strategy to 
pick an appropriate scaling factor and the corresponding 
truncated mode for at least the most commonly used 
types of functions, i.e. the Gaussian type and the super- 
Gaussian type functions. The strategy we are about to give 
only depends on the asymptotic behavior of the function. 
In the scenario where the solution of some differential 
equation needs to be approximated (the exact solution is 
unknown), we could use asymptotical analysis to obtain 
its asymptotic behavior. Thus, our strategy of picking 
the suitable scaling factor is still applicable. A numerical 
experiment is also included to verify the feasibility of our 
strategy. Although it may not be optimal with respect to 
the accuracy, our strategy provides a useful guideline for 
the implementations of HSM. In this paper, the precise 
convergence rate of the HSM to FKE is obtained in 
suitable funciton space by numerical analysis and verified 
by a numerical example. 

Let us draw our attention back to the implementation 
of our algorithm to NLF problems. Through our study 
of HSM to FKE, the off-line data could be well pre- 
pared. However, when synchronizing the off-line data 
with the on-line experiments, to be more specifically, 
updating the initial data according to \\.5) on-line, another 
difficulty arises due to the drifting of the conditional 
density function. The untranslated Hermite functions with 
limited truncation modes could only resolve the function 
well, if it is concentrated in the neighborhood of the 
origin. Let us call this neighborhood as a "window". 
Unfortunately, the density function will probably drift 
out of the current "window". The numerical evidence is 
displayed in Fig. |4.5| To efficiently solve this problem, 
we for the first time introduce the translating factor to 
the Hermite functions and the moving-window technique 
for the on-line implementation. The translating factor 
helps the moving-window technique to be implemented 
more neatly and easily. Essentially speaking, we shift the 
windows back and forth according to the "supporf of the 
density function, by tuning the translating factor. 

This paper is organized as follows. Section II introduces 
the generalized Hermite functions and the guidelines of 
choosing suitable scaling factor to improve the resolution; 
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section III focuses on the analysis of the convergence 
rate of HSM to FKE and a numerical verification is 
displayed. Section IV is devoted to the application to the 
NLF problems. The translating factor and the moving- 
window technique are addressed in detail. Numerical 
simulations of two NLF problems solved by our algorithm 
are illustrated, compared with the particle filter. For the 
readers' convenience, we include the detailed formulation 
of our algorithm in Appendix A and the proof of Theorem 
|2.1| in Appendix B. 

II. Generalized Hermite functions 

Let us introduce the generalized Hermite functions 
and derive some properties inherited from the physical 
Hermite polynomials. 

Let L^(R) be the Lebesgue space, which equips with 
thenorm||-|| = {Jp\-\'^dx)^ and the scalar product {■,■) . 
Let Unix) be the physical Hermite polynomials given by 
Unix) = (-l)"e^'9Je~'^', n e Z and n > 0. The 
three-term recurrence 



■Ho = l, Hi{x) = 2x 
and Hn+iix) = 2x'Hn{x) — 2nHn^iix) 



(2.1) 



is more handy in implementations. One of the well-known 
and useful facts of Hermite polynomials is that they are 
mutually orthogonal with respect to the weight wix) = 

_ 2 

e ^ . We define our generalized Hermite functions as 



h:^^{x) 



=1H„(a(x-/?))e-5"'(-«', (2.2) 



for n £ Z and n > 0, where a > 0, /3 G R are some 
constants, namely the scaling factor and the translating 
factor, respectively. It is readily to derive the following 
properties for \2.2y . 

I) The {H^'^}'^^o forms an orthogonal basis of 
L^{R), i.e. 



/ m'^{x)H'^\x)dx^ 



(2.3) 



where Snm is the Kronecker function. 



2) H"''^{x) is the n eigenfunction of the following 
Strum-Liouville problem 



+ \nU{x) = 0, 



u{x))) 
(2.4) 



with the corresponding eigenvalue A„ — 2a^n. 
3) By convention, H^-^ = 0, for n < 0. Forn G Z 
and n > 0, the three-term recurrence holds: 

2a{x - l3)H^'\x) =^/2^H^'f^[x) (2.5) 
+ .M^HZ'^.ix); 

or 

2a\x - I3)m-^{x) ^^H:l\{x) 

"+l-^n+^l(^)- 



4) The derivative of H"'^ [x) is a linear combination 

of H"J,{x) and H^'^,{x): 



d.HZ'^ix) 

= ^\/A^-H'°fi(3;) - Xn+iH"f^{x) 



aH^L^x)- 



aH^f.ix). (2.6) 



5) Property 1) and 4) yield the "orthogonality" of 

{dxm-f'(x)}^^o- 



dxH"'^ {x)dxH^'^ {x)dx 



4a 



(An + A„+i), if m : 



I = min{n, m}, if \n — m\ = 2; 

0, otherwise. 

(2.7) 

The generalized Hermite functions form a complete or- 
thogonal base in L^(K). That is, any function u G L^(R) 
can be written in the form 



u{x) = ^ u„H"''^{x) 



where {un}^^^ are the Fourier-Hermite coefficients, 
given by 



a 

7^ 



u{x)HZ'^ {x)dx 



(2.8) 



Let us denote the subspace spanned by the first + 1 
generalized Hermite functions by 7?.jv: 

7^iv = span{//o"■''(x),••• , ff°'''(a;)} . (2.9) 

In the sequel, we follow the convection in the asymp- 
totic analysis that a ~ 6 means that there exists some 
constants Ci, C2 > such that C\a < b < C2a; a < b 
means that there exists some constant C3 > such that 
a < Cs,b. Here, C\,C2 and are generic constants 
independent of a, j3 and N. 



A. Orthogonal projection and approximation 

It is readily shown in )2\\ for a > 0, /3 = that the 
difference between an arbitrary function and its orthogonal 
projection onto 7?. at in some suitable function space could 
be precisely estimated in terms of the scaling factor a and 
the truncation mode TV. Let us first introduce the function 
space ;9(K), for any integer r > 0, 



M^a,;3(K) :={iiG L'(R) : \\u\\r,c,.fi < 



fe=0 J 

(2.10) 



where A^ is in l |2.4[ ( and Uk is the Fourier-Hermite 
coefficient in jZ8) . We shall denote W''(R) for short. 
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if no confusion will arise. Also, the norms are denoted 
briefly as j| • jjr - The larger r is, the smaller space W'^(R) 
is, and the smoother the functions in VF'^(R) are. The 
index r can be viewed as the indicator of the regularity 
of the functions. 



Let us define the L^— orthogonal projection P\ 



,a,/3 



L^(R) -i> 7^]v, i.e. given v G -L^(R), 



(2.11) 



The superscript q, /? will be dropped in P^'*^ in the sequel 
if no confusion will arise. More precisely, 

JV 

Pnv{x) := v„H^''^{x), 

n=0 

where Vn are the Fourier-Hermite coefficients defined in 
\2.S\ . And the truncated error j ju— Pjvuj |r, for any integer 
r > 0, has been essentially estimated in Theorem 2.3, 1111 
for a = 1, /3 = 0, and in Theorem 2.1 , 1211 for arbitrary 
Of > and /3 = 0. For arbitrary q > and /3 7^ 0, the 
estimate still holds. 

Theorem 2.1: For any it G VF' (R) and any integer < 
/i < r, we have 



PnuIu < 



■2N- 



(2.12) 



where \u\^ := \ \d^u\\ are the seminorms, if A'^ S> 1. 
The proof is extremely similar as those in 1111 and 1211 . 
Thus, we omit it here and include it in Appendix B for 
the readers' convenience. 



B. Guidelines of the scaling factor 

From Theorem |2.1[ it is known for sure that any 
function in W^'^(R) could be approximated well by the 
generalized Hermite functions, provided the truncation 
mode A'^ is large enough. However, in practice, "suffi- 
ciently" large is impossible. To improve the resolution 
of Hermite functions with reasonable large we need 
the scaling factor a, as pointed out in |4|. Many efforts 
have been made along this direction, refer to (3), l4l . 
1131 , etc. However, the optimal choice of a (with respect 
to the truncation error) is still an open problem. In this 
subsection, we give a practical guideline to choose an 
appropriate scaling factor for the Gaussian type and super- 
Gaussian type functions. 

It is well known that, for smooth functions f(x) = 
'^'^^Q fnH"'^ {x), the exponential decay of /„ with 
respect to n implies that the infinite sum is dominated 
by the first A'^ terms, that is. 



o 



for A'' 3> 1. Thus, the suitable scaling factor is supposed 
to get the Fourier-Hermite coefficients decaying as fast 
as possible. Once the coefficient approaching the machine 
error (say 10~^®), many other factors such as the roundoff 
error will come into play. Hence, it is wise to truncate the 
sequence here. Therefore, we need some guidelines of 



choosing not only the suitable scaling factor a but also 
the corresponding truncation mode A'^. 

Suppose the function f{x) concentrates in the neighbor- 
hood of the origin and behaves asymptotically as e"'''^' 
with some p > and k > 2, as |x| — >■ +00. Our 
guidelines are motivated by the following observations: 

1) The function / decays exponentially fast, as |a;| — > 
00. So, /„ « J_^^ f{x)H"'^ {x)dx, provided L is 
large enough, due to \2.S) . 

2) For the exact Gaussian function e"''^ , p > 0, the 
optimal a is naturally to be \/2p with the truncated 
mode A'^ = 1. It follows by the fact that with 



this choice, e = H^'°{x). thus, e ''^ is 
orthogonal to all the rest of -ff"'", n > 0. That 
is, (e-P^^)o ^ and (e-?^^)^ = 0, n > 1. This 
suggests that the more matching the asymptotical 
behavior of / to e~3" ^ , the faster the Fourier- 
Hermite coefficients decays, and the smaller trun- 
cation mode A' is. 
3) It is natural to adopt the Gaussian-Hermite quadra- 
ture method to compute the Fourier-Hermite coef- 
ficients by \2.9i) . The truncation mode A' has to be 
chosen such that the roots of Hermite polynomial 
Hn+i cover the domain [—aL, aL] where the 
integral \2.S\ is contributed most from both / and 
n = 0,--- ,iV. 
We describe our guidelines for the Gaussian type and 
the super-Gaussian type functions separately as follows. 
Case I. Gaussian type, i.e. f{x) ~ e 

l^l — >■ +00. 

2 122 
1) e-''^ ~ e--" " 



, p > 0, as 
00, which yields 



2) The integrand in i2.8i is approximately e 



Using the machine error 10 to decide the domain 



of interest L, i.e. e 



10" 



yielding that 



L ^ V8p-i In 10; 

3) Determine the truncation mode A' such that the 
roots of Hermite polynomial Hn+i covers approx- 
imately {—aL, aL), where aL « 4Vln 10. 

Case II. Super-Gaussian type, i.e. f{x) ~ e~^^ , as 
^1 — >■ +00 for some fe > 2, p > 0. 

1) Notice that e's"^""^ > e"'"''', when x > 1. Thus, 



we require that e 2" 
that aL ^ \/321nl0; 



10 , which implies 



2) We match e we 2 ^ near x = ±L yields 
that a ^ ySpLt"^ Hence, L x (16p~MnlO)^, 
a fcpfc(lnlO)2 (=; 

3) Determine the truncation mode A' such that the 
roots of Hermite polynomial Hn+i cover approxi- 
mately {—aL,aL). 

To exam the feasibility of our strategy, we explore the 
Gaussian type function f{x) — e~^^ cos(^). Accord- 
ing to the strategy in Case I, we choose the scaling factor 
a ^ ^ 3.1, L ^ y^^^ « 1.9194 and A^ ^ 24. 
As shown in Fig. |2.1[ the truncation error with a = 3.1 
decay the most fast with respect to the truncation mode A'^ 
and approach the machine error near the 20* frequency 
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f(x)=cos(>c;i 0)6"™ , |5=0 



□ □ □ 



□ □ □ □ 



o 

+ 0=3.1 
□ a=1 



10 15 
truncation mode N 



Fig. 2.1. The truncation error v.s. the truncation mode of 

f{x) = cos (^)e-5"=^ is plotted, with ^ = 0, a = 4,3.1 
and 1, respectively. 



mode. Meanwhile, the decay of the truncation error with 
a = 4 and a — 1 are much slower. Moreover, the 
truncation mode A'^ = 24 is appropriate in the sense that 
the next few coefficients start to grow, due to the roundoff 
error. 

Remark 2.1: 1) This strategy is very useful. However, 
it is not the optimal scaling factor a. For example, if 
f{x) — e~2^ , then the optimal scaling factor a — 1 
and A'^ = 0, instead of A'' = 24 from our guideline. 

2) Although the scaling factor helps to resolve the 
function concentrated in the neighborhood of the origin, 
it helps little if the function is peaking away from the 
origin. The numerical evidence could be found in Table 
|4.2| This is the exact reason why we need to introduce 
the translating factor to the generalized Hermite functions 
when applying to the NLF problems, see the discussion 
of translating factor in section IV.B. 

III. Hermite spectral method to ID forward 

KOLMOGOROV equation 
The general ID FKE is in the form 

ut{x,t) =p{x,t)uxx{x,t) + q{x,t)ux{x,t) 

+ r{x,t)u{x,t), for {x,t) eRxR+ 
u{x, 0) = cro(x). 

(3.1) 

The well-posedness of ID FKE has been investigated 
in (2|. We state its key result here. 

Lemma 3.1: (Besala, (2)) Let p{t,x), q{t,x), r{t,x) 
(real valued) together with pa;, p^cx, Qx be locally Holder 
continuous in D — (to,^i) x K. Assume that 

1) p{t, 2;) > A > 0, V (t, x) G D, for some constant 
A; 

2) r{t,x) < 0,y{t,x) G V- 

3) {r-qx+ Pxx){t, x)<0,\/ (t, x) G V. 

Then the Cauchy problem with the initial condition 
u(to, x) = uo{x) has a fundamental solution r{t, x; s, z) 
which satisfies 

< r{t,x;s,z) < c(t- s)~5 



for some constant c and 

/oo 
r(t, x; s, z)dx < 1. 
- 00 

Moreover, if 110(2;) is continuous and bounded, then 

/oo 
T{t, x; to, z)uo{z)dz 
- 00 

is a bounded solution of p.l[ l. 
Through the transformation 

w^x, t) 



□ 



p2 [s,t)ds,t 



(3.2) 



where 

q{x,t) =p~^{x,t) 



<li3:,t) - -p ^Px{x,t) 



+ T 



p ^pt{s,t)d.s 



(3.3) 



equation jSTTJ can be simplified to the following FKE 
with the diffusion rate equals 1 and without the convection 
term. 



Wt{x, t) — Wxx{x, t) + V{x, t)w{x, t), for K X 
'w{x, 0) = ^10(3;), 



(3.4) 



where 

V{x,t) 



1 

+ o / qt{s,t)ds + r[x,t) 



(3.5) 



Remark 3.2: From the computational point of view, the 
form {i.4\ is superior to the original form l |3.1[ ) in general, 
when implementing with the HSM. 

(i) If both the potential V{x,t) and the initial data 



'w{x, 0) are even functions in x, so is the solution to 1 3.4 1 



With the fact that the odd modes of the Fourier-Hermite 
coefficients of the even functions are identically zeros, it 
requires half amount of computations to resolve the even 
functions. 

(ii) Even when V{x, t) and 'w{x, 0) are not even, it is 
still wise to get rid of the convection term, since this term 
will drive the states to left and right, and probably out of 
the current "window". Shifting of the windows frequently 
by the moving-window technique will definitely affect the 
computational efficiency. 



A. Formulation and convergence analysis 

In this subsection, we shall investigate the convergence 
rate of the HSM of solving the FKE. Let us consider 
the FKE j3.4^ with some source term F{x,t). The weak 
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formulation of HSM is to find un{x, t) G TZn such that 

{dtUN{x,t),(p) 

+ {V{x,t)uN{x,t),ip) + {F{x,t),^), 
. un{x, 0) =Pnuo{x), 

(3.6) 

for all if £ TZn- The convergence rate is stated below: 
Theorem 3.2: Assume 

-(i + la;]')^ < V{x,t) < c, 

for all {x,t) G K X (0,r), for some 7 > and 
some constant C. If uo G VK'^(K) and it is the so- 
lution to 1 3.4 1 with source term F(x,t), then for u G 
L°°{0,T;W''(R)) n i^(0,r;'tV(K)) with r > 27 and 

r 47-21- + 2 . 1 

iV>max<^ci 2^-1 max{(Q/3) 



Proof of Theorem \3.2\ Denote Un ~ Pnu for 



a -T max{(a^)*^,l} 



it yields that 



\\u - UN fit) < ca-^''-^ ma.yi{(apf'' , l}N^''-\ 

(3.7) 

where c* depends only on T, |!itjjioo(o,T;W"-(E)) and 

11^111-2(0, T;VK '-(«))• 

Before we prove Theorem |3.2[ we need some estimate 
on ||x'^i9J!2j((^j;)||2^ fQj. jijjy integer ri,r2 > 0: 

Lemma 3.2: For any function u G W^^^^^ (M.), with 
some integer ri , r2 > 0, we have 



\x''^ d^'^ u\f < a ^'■^ ^ max{(Q^)^''\ 



(3.8) 



Proof: For any integer ri , r2 > 0, 



r-i I I 2 

X Oj; lt|| = 



y^^u„x''^d^'^Hn''^{x) 

^ 00 r2+ri 



by ([23} and ( |2.6[ l, where for each n fixed, a„,fe is a 
product of 2(ri + r2) factors of or ^A„+j, with 
— r2 — ri < j < r2 + ri. Let n* > such that 



a /? ~ And notice that A„+_, ~ A^+i for 

, T- 1 O Tr\r n _L 



n + j > and W^fAx) = for n + j < 0. Hence, 



we have 



00 



Jl — Tl * + 1 

-2ri-l o\2ri 



<Q max{(a/3) ''^ , 1}| juj 1^^+^ 



simplicity. By i3.4i with source term F{x,t) and the 



definition of Un, we obtain that 

= {dt{u - Un),<p) = - {ux,'fix} + {V{x,t)u,ifi) 
+ {F{x,t),(f) - {dtUN,'p) 

^ {dtUN,'p) = - {Ux,'Px} + {V{x,t)u,if) 

+ {F{x,t),^), (3.9) 
for all (/p G 72. AT. Combine with l |3.6| l, it yields that 

{dt{uN - Un), f) ^ - (dxiuN - u), ^x) 

+ {V{x,t){uN — u),ip), 

for all </? G 7?.]v. Set qn = un — Un- Choose the function 
ip — 2qn, then we have 

dtWoN^ = - 2\\dxQNf - 2{dx(UN - u), BxQn) 
+ 2{V{x,t)QN, Qn) 

+ 2{V{x,t)[UN ~u),qn). (3.10) 
By Young's inequality, 

\{dx{UN ~ u),dxQN)\ 

< l\\dx{UN -u)f + \\dxgNf. 

(3.11) 

The assumption V{x,t) < C for {x,t) G M x (0,T) 
yields that 

{V{x,t)gN,gN) <C\\eNf, (3.12) 
for (x,t) G R X (0,T). Moreover, we have 

\{V{x,t){UN-u),QN)\ 

< ^\\V(UN^n)f + l\\gNf, 
(3.13) 



by Cauchy-Schwartz's inequality. Substitute l |3.11[ l-( |J. 13[ l 
into l |3.10| l, we obtain that 

dt\\gN\\' -{C+l)\\gNf 

<\\V{Un -u)f + ^\\dx{UN - n)f. (3.14) 

Notice that V > — (l+|a;p)^, for some 7 > 0. Essentially 
by the estimate in Lemma [^2] we can estimate 

\\V{UN-U)f 

<||(1 + \x\'')-'{UN-u)f < Wix^-' + 1){Un - u)\\'' 

00 

<a"'*^"^max{(a^)''^,l} ^ ^t'iiii-l 

n = JV + l 

+ \\Un — uf 
<a""-'-' miix{iaP)*\l}N^^-''\\u\\l 



I — 2r — Iat — '"ll l|2 

+ a N kt r- 



(3.15) 



for any integer r\,r2 > 0, by l|2.3 



The estimate of the second term on the right hand side of 
l|3.15^ is followed by Theorem |2.1| Again, by Theorem 
|2.1| we obtain 

\\dx{UN - u)f = \Un - u\i < a-^'-+'N'-'-\\u\\l 

(3.16) 
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Substitute ([STTSj, \3.16\ into l |3.14^ , we obtain 

dtWQNlf -iC+l)\\gN\\'' 

provided that 

Af>max<^ci ^^-i max {(q/3) 1} ^-^^ , 
a^-T max{(a/3)*^,l}"5^ 



/' 

^0 



Therefore, it yields that 

\\QN\\\t) <a-^^-Vax{(a/3)^^l}iV2^-'' 

'0 

By the triangular inequality and Theorem |2.1| 

||u — Mjv||^(t) 

< ||£)Ar||' + llM-[/ivl|' 

<a-^->-^N^'^-^ [Ml 

+ max{(a/?)'''',l} /* e-<^+i'<*-"'||u||?,(s)ds 

< c'a-"^'^ max{(a/?)*^, IjA^^^"", 

where c* is a constant depending on | |loo(o,t;W(m))> 



11,2(0, T;W"-(E)) 



and T. 



S. Numerical verification of the convergence rate 

To verify the convergence rate of HSM, we explore an 
ID FKE with some source F{x,t). The exact solution 
could be found explicitly as our benchmark. The error 
v.s. the truncation mode is plotted. 

We consider the ID FKE 



ut = u-xx — X u + (sin t + cos t + 3x)e 



u{x, 0) = xe 



(3.17) 



for (a;, t) £Rx [0, T]. It is easy to verify that u{x, t) = 
[x + sint)e~ 2^ is the exact solution. 

Notice that the initial data, the potential and the source 
in l |3.17^ are all concentrated around the origin. So, we 
set /3 = 0. For notational convenience, we drop /3 in this 
example. As to the suitable scaling factor q, from our 
strategy in section II.B, we know that it is better to let 
a = 1. However, if we do so, the first two modes will 
give us extremely good approximation. Hence, the error 
v.s. the truncation mode won't be perceivable. Due to this 
consideration, we pick a = 1.4 (a little bit away from 1, 
but not too far away so that it won't affect the resolution 
too much). The formulation l|3.6^ yields 



{dtUN,(p) = — {dxUN,dx(p) — {xUN,X(p) 

+ {F{x,t),^), (3.18) 

for all ((9 G TZn- Take the test functions tp = H"(x), n = 
0, 1, • • • , A'^, in l|3.18|l. The numerical solution ujv £ TZn 




20 25 30 35 

truncation mode N 



Fig. 3.2. The L^.g^-ors of the HSM to FKE pi?) v.s. the 
truncation mode N = 5, 15, 25, 35 and 45 is plotted, with a = 
1.4, /3 = and the time step dt = 10~^. 



can be written in the form 



UN 



The matrix form of ( |3.18^ follows from l |2.5[ > and \2J\ : 

dta{t) = Aa{t) + f{t), (3.19) 

where a(t) = (ao(t), ai(t), ■ • ■ , ajv(t))^, /(t) = 

^/o(t), fi{t), ■ • ■ , ,/]v(t)^ are column vectors with A'^ + 

1 entries, fi{t), i = 0, 1, • • • , A^, are the Fourier-Hermite 
coefficients and ^4 is a penta-diagonal (A'' + 1) x (A'^ + 1) 
constant matrix, where A = —Ai — A2, 

2 



fc = min{j,j}, |i-j|=2, 



, 0, otherwise, 



and 



{ V(fc + l)(fc + 2) 
2q2 

fc = min{j,j}, \i-j\=2, 

(2i + 1) . . 

* =3, 



2q2 
I, otherwise. 

The errors v.s. the truncation mode A'^ at time T = 
0.1 is plotted in Fig. |3.2| The ODE l|3TT9j is numerically 
solved by central difference scheme in time with the time 
step dt = 10"^. It indeed shows the spectral accuracy of 
HSM. 

IV. Application to nonlinear filtering 

PROBLEMS 

Recall the brief description of our algorithm in the 
introduction (and more details in Appendix A), the off- 



7 



HSM to ID FKE and its application to the NLF problems 



line computation is to numerically solve the FKE ( |1.4| l 
repeatedly on each interval [ri,ri+i]. Equation l |1.4| is in 
the form of ( |3.1| l with 



r(x,t) 



\h^/S + Q{gl +ggxx) - fx 



where Q, S, /, g and h are in {OJ. 

A. Existence and Uniqueness of the solution to the FKE 

In this subsection, we interpret the well-posedness 
theorem (Lemma 1 3. 1^ for general ID FKE in section III 
in the framework of the NLF problems. 

Proposition 4.1 (Existence): Let /, g, h in are 
Holder continuous functions m T) := R x {ta,ti). gx, 
gxx and fx exist and are also Holder continuous in T>. 
Assume further that 



1) 



> A > 0, for some A > 0; 



2) S > 0; 

3) -\h^/S - fx + Q{gl + ggxx) < C, for some 
constant C, 

for (x, t) G T). Then there exists a bounded solution 
u{x, t) to l |3.1^ , if the initial condition uo{x) is continuous 
and bounded. 

Proof: Conditions l)-3) in Lemma [3TT] are directly 
translated into conditions l)-3) in this proposition with 
C < 0. For C > 0, let v{x,t) = e-'^'*~*»)'u(a;, f), then 
V satisfies 

Vt (a;, t) =p{x, t)v^^ {x,t)+q{x, t)v^ {x, t) 

+ {r{x,t)~C)v{x,t), (4.20) 

for {x,t) G T), with the initial condition v{x,to) = 
110(2;). The coefficients of l |4.20^ satisfy the conditions in 
Lemma [3T| Thus, we apply Lemma [3TT] directly to l |4.20^ . 
The existence of the solution to \2>A\ follows immediately. 

■ 

Remark 4.3: In practice, the initial data of the condi- 
tional density function is either compactly supported or 
exponentially decays as \x\ — >■ +00. So, the assumption 
on the initial data in Proposition |4.1| is valid. 

For concise of notations, we give the uniqueness for 
the equation instead of l |3.1| >. It can be easily trans- 
formed into each other, due to the bijective transformation 

Proposition 4.2 (Uniqueness): There exists 
a unique solution to l |3.4| ( in the class that 
{u : lim\x\-foo uux = 0} if V{x,t) is bounded 
from above in V. 

Proof: Case I: Assume V{x,t) < in 2?. Suppose 
there exist two distinct solutions to l |3.4^ , say ui and U2. 
Denote rj :— ui — U2, and 77 satisfies 



Vt = Vxx + V{x,t)r], 



(4.21) 



in 2? with the initial condition rj{x,to) = 0. Use the 
standard energy estimate, i.e. multiplying ( |4.21^ with rj 
and integrating with respect to a; in R: 

1, 



/ 



T^xW + / V(x,t)r]''dx < -\\Vx\\ < 0, 



TABLE 4.2 

Trunction error v.s. the "peaking" po of the 
Gaussian function f{x) 



PO 


erroro 


errors 


-1 


3.3 X 10"" 


1.1 X IQ-^ 





8.2 X 10-1'' 


7.7 X IQ-*' 


1 


1.6 X 10"^^ 


1.8 X 10"^ 


2 


1.8 X IQ-^ 


3.3 X IQ-'-'-^ 


3 


7.7 X IQ-" 


8.2 X 10"-^^ 


4 


1.1 X IQ-^ 


1.6 X IQ-'-'-^ 



The scaling factor is chosen to be 1 according to the guideline in 
section II. B and the truncation error is N = 24. The truncation 
en'ors with different translating factor /3 is denoted as error^, 
which is defined as ||/ - Yln^o UH^'^^W- 



by the 
lim 



that 



integration by parts, and the facts that 
> rjrjx = and V{x,t) < in D. This yields 



ll'7ll'W< 



for t £ {to,ti). With the fact that r){x,tQ) = 0, we 
conclude that 77 = in D, i.e. ui =112. 

Case II: Assume V{x,t) < C, for some C > 0. 
We use the strategy in the proof of Proposition |4.1[ Let 



v{x,t) = e~^^ °'u{x,t), then v satisfies i3.4i with the 
potential V{x, t) — C < in V. By case I, we conclude 
the uniqueness of v, so does u. ■ 
Remark 4.4: The similar conditions as in Proposition 
|4. 1| are derived to guarantee the well-posedness of the 
"pathwise-robusf DMZ equation in |22| and to estab- 
lish the convergence of our algorithm in II7I . They 
essentially require that h has to grow relatively faster 
then /. They are not restrictive in the sense that most 
of the polynomial sensors are included. For example, 
f(x) = fox\ g{x) = go(l + a;^)''' and h{x) = hox'- , 
with 5, Q > 0, /o,<7o and fto are constants, j,k,l G N, 
provided I > max { , 2fc — 1 } . 



B. Translating factor P and moving-window technique 

As we mentioned in the introduction, the untranslated 
Hermite functions with suitable scaling factor could re- 
solve functions concentrated in the neighborhood of the 
origin accurately and effectively. However, the states of 
the NLF problems could be driven to left and right during 
the on-line experiments. It is not hard to imagine that the 
"peaking" area of the density function escapes from the 
current "window". 

The translating factor /3 is introduced under the cir- 
cumstance that the function is peaking far away from the 
"window" covered by the current Hermite functions. We 
translate the current Hermite functions to the "support" of 
the function, by letting the translating factor /3 near the 
"peaking" area of the function. 

In Table 14.21 we list the truncation error of the Gaus- 



sian function f{x) = e 2 



with various po = 



— 1, ■ ■ ■ ,4 and different translating factors /? = or 3. 
According to the guidelines in section II.B, the scaling 
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u/=o i^L^, + I3j,l3j + L„). For simplicity and clarity, 
let us assume further that 

1) The operator (i — ^h'^S~^h) is not explicitly 
time-dependent; 

2) The time steps are the same, i.e. r^+i — Ti — At. 
For the storage of the off-line data, on each interval 
{—L^+Pjjfij+Lw), it requires to store (TV + 1)^ 
floating point numbers. Hence, the total (J + 1) intervals 
requires to store (J + 1){N -\- Vf' floating point numbers. 
As to the number of the flops in the on-line computations, 
if no moving-window technique is adopted during the 
experiment, for each time step, it requires 0{{N + 1)^) 
flops. The number of the flops to complete the experiment 
during [0,T] = u'^-^[Ti,Ti+i] is 0{k{N + lf). Suppose 
the number of shifting the windows during [0, T] is P, 
then the total number of flops is O {{k + P){N + 1)^). 

Remark 4.5: If either assumption 1) or 2) is not satis- 
fied, then the real time manner won't be affected. That is, 
the number of the flops in the on-line experiment remains 
the same. But the off-line data will take more storage 
as the trade-off. To be more specific, on each interval 
(— L„ + /3j, Pj + Lw), it requires to store fc x (A'' + 1)^ 
floating point numbers, where k is the total number of 
time steps. Therefore, the total storage is k{J+l){N+l)'^ 
floating point numbers. 

C. Numerical simulations 

In this subsection, we shall solve two NLF problems 
by our algorithm: the almost linear sensor and the cubic 
sensor. Since the drift term could always be absorbed 
into the potential V{x) by the transformation \3.2\ , for 
simplicity, in our examples, we set / = 0. Our algorithm 
is compared with the particle filter in both examples. 
The particle filter is implemented based on the algorithm 
described in 1 1 1. And the systematic resampling is adopted 
if the effective sample size drops below 50% of the 
total number of particles. As we shall see, our algorithm 
surpasses the particle filter in the real time manner. 

1) Almost linear filter: We start from the signal obser- 
vation model 

dxt = dvt 

{Pj}j=o, for some J > 0, need to be pT^ared before- I '^V* = + 0.25co8Xt)dt + dwt, 

hand, such that [—L,L] C u/=o {—L^ + j3j,l3j where Xt, yt € R, vt, wt are scalar Brownian motion 



factor is a = 1 and the truncation mode = 24. As 
shown in the table, the further the function is peaking 
away from the origin, the larger the error is with untrans- 
lated Hermite functions. But with appropriate translating 
factor, the function could be resolved very well with the 
same scaling factor, for example, error,-? ~ 10"" for 
f{x) =e-3{— 3)'. 

Indeed this fact motivates the idea of moving-window 
technique. The suitable width of the window could be pre- 
determined if the trunction error of the density function 
v.s. various "peaking" po is investigated beforehand. To be 
more precise, suppose we know the asymptotical behavior 
of the density function of the NLF problem from the 
asymptotical analysis, say ~ e"''^ , with some p > 0, 
k > 2. According to the guideline in section II.B, the 
suitable scaling factor a and the truncation mode with 
/3 = could be chosen. With these parameters, the similar 
table as Table |4.2| could be obtained, i.e. the truncation 
error (erroro) of the function e~*''^~''°' v.s. various po - If 
the error tolenrance is given, then the appropriate width of 
the window is obtained according to the table. Take Table 
|4.2| as a concrete example. If the asymptotical behavior 
of the density function is e~2^ , then the scaling factor 
a — 1 and the truncation mode — 24. Suppose we 
set the error tolerance to be 10~^, then the suitable width 
of the window would be 3 + 3 = 6, from the first two 
column of Table [4. 2 1 The window covers the origin would 
be [-3,3]. 

Our algorithm with moving-window technique is il- 
lustrated in the flowchart Fig. |2.6| It reads as follows. 
Without loss of generality, assume that the expectation 
of the initial distribution of the state is near 0. During 
the experimental time, say [0,r], the state remains in- 
side some bounded interval [— L,L], for some L > 0. 
We first cover the neighborhood of by the untrans- 
lated Hermite functions {H"'^^}^^q, where a, N can 
be chosen according to the guidelines in section II.B. 
With the given error tolerance, the suitable width of 
the window could be pre-defined, denoted as L„. If 
[—L,L] C [— Lu,,Lu,], then no moving- window tech- 
nique is needed. Hence, the on-line experiment runs 
always within the left half loop in Fig. |2.6| Otherwise, 



The off-line data corresponding to different intervals 
{—Lw+I3j,l3j+Lw) have to be pre-computed and 
stored ahead of time. During the on-line experiment, if the 
expectation of the state E[xt] moves accross the boundary 
of the current "window" (the condition in the rhombic box 
in Fig. |2.6| is satisfied), the current "window" is shifted 
to the nearby window where E[xt] falls into. That is, the 
right half loop in Fig. |2.6| is performed once. 

Let us analyze the computational cost of our algorithm. 
Notice that only the storage capacity of the off-line data 
and the number of the flops for on-line performance 
need to be taken into consideration in our algorithm. 
Without loss of generality, let us assume as before 
E[a;](0) is near and our state is inside [~L,L] C 



processes with E[dVf dvt 



1, E[dwt dwt 



Suppose the signal at the beginning is somewhere near 
the origin. 

The corresponding FKE in this case is 



Ut = 



1 _1 2 

^u^. ^x 



(1 + 



COS a; u 



(4.22) 



Assume further that the initial distribution of xq is 
uq{x) — e 2 . This assumption is not crucial at all. The 

non-Gaussian ones, for example ito(a;) — e 2 , vvill give 
the similar results as the Gaussian one. 

It is easv to se e that the asymptotical behavior of the 

solution to 1 4.22 1 is e 2 . With the guidelines in section 
II.B, we choose a = 1, /3 = and TV = 25 for the 
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Almost linear sensor 



14 




-2 • • • • • • • • • 1 

5 10 15 20 25 30 35 40 45 50 
time 

Fig. 4.3. Almost linear filter is investigated with our algorithm 
and the particle filter with 10 and 50 particles. The total experi- 
mental time is T = 50s. And the update time is At = 0.01. 



starting interval. We shall run the experiment for the total 
time T — 50s. Thus, we expect the density function 
probably will move out of the starting interval. Table 
|4.2| suggests that the appropriate width of the window 
should be 3, if the error tolerance is set to be lO^"". 
We shall overlap the adjacent windows a little bit to 
prevent frequent shifting of windows. Let us take the 
width of the overlaped region to be 0.5. Therefore, as the 
preparation for the moving-window technique, we shall 
prepare the off-line data for [-19.5,-13.5], [-14,-8], 
[-8.5,-2.5], [-3,3], [2.5,8.5], [8,14] and [13.5,19.5]. 
The correpsonding /3s are -16.5,-11,-5.5,0,5.5,11 
and 16.5. The barrier in the rhombic box in the flowchart 
Fig |2.6| should be 3 (the width of the "window"). 

Our algorithm is compared with the particle filter with 
10 or 50 particles in Fig. |4.3| for the total experimental 
time T = 50s. The time step is At = 0.01s. All three 
filters show acceptable experimental results. It is clear 
(between time 10 to 30) that the particle filter with 50 
particles gives closer estimation to our algorithm than that 
with 10 particles. But as to the efficiency, our algorithm 
is superior to the particle filter, since the CPU times of 
particle filter with 10 and 50 particles are 5.00s and 
35.75s respectively, while that of our algorithm is only 
2.62s. As to the storage, the size of the binary file to keep 
the off-line data is only 35.5kB. During this particular on- 
line experiment, the window has been shifted for 13 times, 
which can't be perceivable from the figure at all. And it 
seems that the moving-window technique doesn't affect 
the efficiency of our algorithm. 

2) Cubic sensor in the channel: We consider cubic 
sensor in the channel xt G [—3,3]: 

dxt ~ dvt 

3 (4.23) 
dyt — Xfdt + dwt, 

where xt, yt G K, vt, Wt are scalar Brownian motion 
processes with Eldvfdvt] = 1, Elduufdiut] — 1. 
Assume the initial state is somewhere near 0. 



Cubic sensor 




5 10 15 20 25 30 35 40 45 50 
time 



Fig. 4.4. Cubic sensor in the channel is experimented for T = 
50, with the time step At = 0.01s, by both particle filter and 
our algorithm. 




50 



X 



Fig. 4.5. The normalized density functions are plotted every 
other Is for the cubic sensor in the channel. 

The FKE is 

1 1 6 

Ut = -Uxj: — -X u. (4.24) 

Furthermore, we assume the initial distribution is uq [x) = 
g-i /4 gjjj^-g {jjg gf^fg inside the channel, we set 
our translating factor /3 = and the moving-window 
technique won't be used. According to section II.B, we 
choose the scaling factor a ^ 23 (iB^) ' ^ 2.4637, and 
the truncated mode TV ~ 45. 

In Fig. |4.4| we compare our algorithm with the particle 
filter with 50 particles for T — 50s. The observation data 
come in every 0.01s. Fig. |4.4| reads that both filters work 
very well. The result of our algorithm nearly overlaps with 
that of the particle filter all the time. However, the CPU 
time of our algorithm is 4.90s, while that of particle filter 
is 37.17s. With our algorithm, the on-line computational 
time for every estimation of the state is around 0.001s, 
which is 10 times less than the update time 0.01s. This 
indicates that our algorithm is indeed a real-time solver. 
The normalized density functions, which is defined as 
"^^''l — TT, have been plotted every other Is in Fig 

max^gj u(x.t] f J o 

El 
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V. Conclusions 

In this paper, we first investigate the HSM applied to the 
ID FKE. It is well-known that the choice of the scaling 
factor Q is crucial to the resolution of HSM. We give a 
practical guidelines to help choosing the suitable one. The 
convergence rate of the HSM has been shown rigorously 
and has been verified by a numerical experiment. As an 
important application, we solve the NLF problem, by 
using the algorithm in fTTl . in the last section, where 
solving ID FKE serves as the off-line computation. To 
capture the state even if it drifts out of the "window", 
translating factor of Hermite functions and the moving- 
window technique are introduced. The translating factors 
help the switch of the windows back and forth easier, 
according to the "support" of the density function of the 
state. We analyzed the computational complexity of our 
algorithm in detail, with respect to the storage capacity 
of off-line data and the number of flops of the on-line 
computations. Finally, two online experiments - almost 
linear filtering and cubic sensor in the channel - are 
reported. The feasibility and efficiency of our algorithm 
are verified numerically, which surpasses the particle filter 
as a real time solver. 

Appendix A 
The detailed formulation of our algorithm 

Starting from the signal model ( |l.l[ l, the DMZ equation 
is derived for the unnormalized density function 
a{x, t) of the states Xt conditioned on the observation 
history Yt = {j/s : < s < f}. In real applications, one 
is more interested in the robust state estimators. Hence, 
for given observation path yt, let us make an invertible 
exponential transformation 

a{x,t) ^ explh^ {x,t)S~'^ {t)yt]p{x,t). (A.l) 

The "pathwise-robust" DMZ equation is obtained: 

f^ix,t) + ^^{h^s-YytPix,t) 



the robust DMZ equation with yt 

Ti-i < t < Ti, i — 1,2, - ■ ■ ,k 



on the interval 



= exp(-/i^S ^yt) 

■ [exp {h^ S''^yt)p{x,t)] 
, pix,0) ^ ao{x). 



2 



(A.2) 



The exact solution to l |A.2| l, generally speaking, doesn't 
have a closed form. Hence, we developed an efficient 
algorithm to construct a good approximation in 1 17|. 

Let us assume that we know the observation time 
sequence = ro < ri < ■ • ■ < Tk = T apriorily. 
But the observation data {yr;} at each sampling time r,;, 
i = 0, ■ ■ ■ ,k are unknown until the on-line experiment 
runs. We call the computation "off-line", if it can be per- 
formed without any on-line experimental data (or say pre- 
computed); otherwise, it is called "on-line" computations. 
One only concerns the computational complexity of the 
on-line computations, since this hinges the success of "real 
time" application. 

Denote the observation time sequence sls Vk ~ {0 = 
To < Ti < ■ ■ ■ < Tk = T}. Let Pi be the solution of 



2 



^{x,t) + ^^{h^S-^y yr,_,p,{x,t) 

= exp (^-h^ S~^yT,_-^ 

■ \eyip(h^ S~^yr^_^^^ p,{x,t)^ 
pi(a;,0) = o-o(a;), 

or 

_ Pi{x,Ti-i) = pi_i(a;,Ti_i), for i = 2, 3, • 



Define the norm of Vk by \'Pk\ — supi<;<j.(ri 
Intuitively, as \'Pk\ — > 0, we have 

k 



in some sense, for all < t < T, where p{x,t) is the 
exact solution of ( |A.2| (. To maintain the real time manner, 
our algorithm resorts to the following proposition. 

Proposition 1.3: For each r^-i < t < Ti, i = 



■ ,k. 

(A3) 

Ti-l). 



,■■ ■ ,k, pi{x, t) satisfies i A.3 i if and only if 

,{x,t) = exp \h^ {x,t)S'^{t)yr^_,'^p,{x,t), (A.4) 



satisfies the FKE ( fL4| l. 

The initial data need to be updated as \i.5\ , followed from 
JA31. 



With the observation time sequence known {ri}*Li, 
we obtain a sequence of two-parameter semigroup 
{U{t,Ti-i)}*^^i, for Ti-i < t < Ti, generated by 
the family of operators {L — |/i"^5'^^/i}t>o. The off- 
line computation in our algorithm is to pre-compute 
the solutions of \\A) at time t = r^+i, denoted as 
{U{Ti+i,T,)(l)i}fl^, where {<;/>; (a::)}gi (chosen as the 
initial data at i = Ti) is a set of complete orthonormal base 
in L^(R"). These data should be stored in preparation of 
the on-line computations. 

The on-line computation in our algorithm is consisted 
of two parts at each time step r^-i, i = 1, • • • ,k. 

• Project the initial condition Ui{x,Ti-i) G L^(R") 
at t = Ti-\ onto the base {(j)i(x)}^-^, i.e., 
Ui{x,Ti-i) = '^^^Ui^i4>i{x). Hence, the solution 
to (|1.4|l at t = Ti can be expressed as 



(A.5) 



Ui{x,Ti) ^U{Ti,Ti-l)Ui{x,Ti-l) 

OO 

= [U{Ti,Ti-l)(jll{x) 



where {U{Ti,Ti^i)(l}i{x)}fZi have already been 
computed off-line. 

Update the initial condition of ^IA\ at n with the 
new observation y^. Let us specify the observation 
updates (the initial condition of ^IA\ ) for each 
time step. For < t < ri, the initial condition 
is ■ui(x,0) = a-o{x). At time t = ti, when the 
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observation is available, 

U2{x, Ti) 

^ exp [h'^(x, Tl)S~^{Tl)yrj^]p2{x, n) 

exp [h''^{x,Ti)S^^{Ti)yri]ui{x, n), 

with the fact yo = 0. Here, ui{x,ti) = 
[W(Ti,O)0i(a;)], where {uij}fZi is com- 
puted in the previous step, and {U{Ti,0)(j)i{x)}'j^i 
are prepared by off-line computations. Hence, we 
obtain the initial condition U2{x,ti) of for the 
next time interval ri < t < T2. Recursively, the 
initial condition of l|1.4b for r,;_i < t < is 



and 



Ui{x,Ti^l) =ex.p[h'^ {x,Ti-l)S ^(Ti_l) 

(yr,_i - 2/r,_2)]Wi-l(2:,n-l), 

(A.6) 

for i — 2,3, ■•■ where Ui_i (a;, Ti_i) = 

Ui-2,l pl{Ti^l,Ti^2)<t)l{x)]. 

The approximation of p[x,t), denoted as p{x,t), is ob- 
tained 

k 

P{^,t) = ^XlTi_i,Ti]it)pt{x,t), (A.7) 

i = l 

where pi{x,t) is obtained from Ui{x,t) by l|A.4[. And 



a{x,t) could be recovered by i A.l 



Appendix B 
The proof of Theorem I2.1I 

Proof of Theorem \2.1\ By induction, we first show 
that for p = Q. For any integer r > 0, 

y OC 

II D l|2 V'^ .2 

Q ^ — ' 

n = JV + l 



y OC 

yT^ >-r -.r ~2 



n = N + l 



(B.l) 



Suppose for 1 < < r, ( |2.12^ holds for ^ — 1. We need 
to show that j2.12| ( is also valid for p. It is clear that 

\U — Pnu\^ <\dxU — PnOxU]^-! 

+ \PNdxU - dxPNu\^^x. (B.2) 

On the one hand, due to the assumption for /i — 1, we 
apply j2.12[ l to dxU and replace p and r with p ~ 1 and 
r — 1, respectively: 

\dxu - PNdxu\^-i <a''-''-^N'^\\dxu\\r-i 

<a""'-5iV^||ti||^, (B.3) 

where the last inequality holds with the observation that 

°° 2 



a 



: / udxH^" 
Jm 



20F 
2^ 



(a;)cia; 

uHl'J^^{x)dx, by (12.6 1 



n A 

1tn + l 7, — fin-l- 



2 ^ 2 

On the other hand, by the virtue of \2.()\ 

PnQxU — OxPnu 

oo 

^PN'^UndxH^''^ {x) - UndxH"''^ {x) 

n=0 n=0 

n=0 
^ iV + 1 

n = 

-2 5Z \/An + lUnH"+i(s) 
n=0 

1 ^ 

+ 2 

= ^\/Ajv+i |^'Ujv-H'^'+i(3;) + ujv+i//^'''(x)j 
This yields that 

IPnOxU — dxPNu\l,-l 

<Xm+i {ul\H^iA^)\l.i+ul+,\H%-^{x)\ 



(B.4) 



due to the property of seminorms. Moreover, we estimate 

ul and \H^-'^{x)\l_-^, for k = N,N + l: 

OO 

ul! < ill < - Pn-iu\\^ < a~^''N'''\\u\\l, 

(B.5) 



by ||rT}. Similarly, -u^^+i < Q-^'^iy'llull^. And 

<a'^\\H'^''^{x)\\l_i, by Lemma [T2| 

— l-i/i — 1 ^- — — 1 

=Q! Ajv < a Ajv+i, 
since {H'^'^)^ = 5fejv, for fc £ 2,'^. Simi larly, 



(B.6) 



-H'^+iIm-i ^ «"^A^ ,\. Substitute (|B.5l and (|B.6|( into 



B.4^ , we get 

IPnOxU - dxPNu\l-l <a~^''~^A*'"''A'i^_|_il|M||r 

<^2M-2r-l^M-r||^||2^ (B.7) 



by the fact that \n = 2Na^. Combine jB.2| , |B3 
(|B.7^, we arrive the conclusion. 



and 
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Start 



Compute the Her- 
mite coefficient 

{Mi,n(Tj)}^=0 

the base {H-'^{x)}^^^ 



Update, i.e. evaluate 
Uiix,Ti) at +/3}^^o 




Input Q, 0, and 
initial data uo 



Compute E[,x]{0) 
and 




Evaluate Ui{x,Ti+i) 




Evaluate Ui(x, Tj) 

using {Ui,niTi)}^^o 



Compute Hermite 
coefficient {ui,n{Ti)}'^^o 
under the shifted 

bases _^ 



Let/3 = 



Synchronized with 
the corresponding 
part of the off-line 
data, i.e. the values of 

{U(Tt+un)H^-P(x)}';^^^ 



Fig. 2.6. The flowchart of our algorithm, where /3' e {0j}j=o- 



14 



